Reading, processing, and exploring Visium data

Objectives

• Quality control and post-processing of a Visium dataset
• Learn visualization approaching for spatial data
• Apply descriptive statistics for lattice-based spatial data

---

This notebook will analyze an anaplastic Wilms Tumor sample from a 3 year old male biopsied at initial diagnosis. The dataset comes from the Single-cell Pediatric Cancer Atlas Portal, project SCPCP000006.

Set up

Load libraries and set random seed:

suppressPackageStartupMessages({
  library(SpatialExperiment) # core object 
  library(VisiumIO) # I/O, but we'll still use :: in the code
  library(ggspavis) # plotting
  library(ggplot2) # plotting
  library(patchwork) # plotting
  library(Voyager) # SFE object & spatial descriptive stats
})
set.seed(2026)

Paths:

data_path <- here::here( 
  "data", 
  "SCPCS000190",
  "SCPCL000429_spatial"
) 

What’s in Space Ranger output? Must note that this was processed with an older version of Space Ranger 1.3.1.

dir(data_path)

[1] "filtered_feature_bc_matrix" "raw_feature_bc_matrix"     
[3] "spatial"                   

filtered_feature_bc_matrix
raw_feature_bc_matrix
spatial

dir(file.path(data_path, "spatial"))

[1] "aligned_fiducials.jpg"     "detected_tissue_image.jpg"
[3] "scalefactors_json.json"    "tissue_hires_image.png"   
[5] "tissue_lowres_image.png"   "tissue_positions_list.csv"

aligned_fiducials.jpg
detected_tissue_image.jpg
scalefactors_json.json
tissue_hires_image.png
tissue_lowres_image.png
tissue_positions_list.csv

Import Visium data

We’ll use the VisiumIO package:

# define handle, essentially
tenx_object_raw_matrix <- VisiumIO::TENxVisium(
  resources = file.path(data_path, "raw_feature_bc_matrix"),
  spatialResource = file.path(data_path, "spatial"),
  # this argument is needed: the default argument includes "cytassist" 
  # but ScPCA samples don't have that one; only have these
  images = c("lowres", "hires", "detected", "aligned")
)

# actually import it
spe <- VisiumIO::import(tenx_object_raw_matrix)

Let’s have a look:

spe

class: SpatialExperiment 
dim: 60664 4992 
metadata(2): resources spatialList
assays(1): counts
rownames(60664): ENSG00000284662 ENSG00000186827 ... ENSG00000277475
  ENSG00000275405
rowData names(3): ID Symbol Type
colnames(4992): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
  TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
colData names(4): in_tissue array_row array_col sample_id
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):
spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
imgData names(4): sample_id image_id data scaleFactor

It’s essentially an SCE, but with a few more bells and whistles.

Note that we have 4992 spots - this is not a random value! It’s the number of spots on Visium’s 6.5 x 6.5 slide.

Unlike SCEs, we have spatial information in a dedicated spatialCoords slot:

# handy function from our code package SpatialExperiment
# literally is giving x/y coordinates
head(spatialCoords(spe))

                   pxl_col_in_fullres pxl_row_in_fullres
AAACAACGAATAGTTC-1               5084               2772
AAACAAGTATCTCCCA-1              14913              12629
AAACAATCTACTAGCA-1               8159               3352
AAACACCAATAACTGC-1               5476              14452
AAACAGAGCGACTCCT-1              13971               5505
AAACAGCTTTCAGAAG-1               4324              11289

We also have a slot that just holds the images, imgData, although looking at it directly isn’t very interesting. But, this stored information will help us make figures!

imgData(spe)

DataFrame with 4 rows and 4 columns
    sample_id    image_id   data scaleFactor
  <character> <character> <list>   <numeric>
1    sample01      lowres   ####   0.0286028
2    sample01       hires   ####   0.0953425
3    sample01    detected   ####          NA
4    sample01     aligned   ####          NA

We should also point out what’s in colData. Unlike an SCE, this isn’t cell metadata. In an SPE, it’s spot metadata; recall that each spot may contain multiple cells, which can even include partial cells!

head(colData(spe))

DataFrame with 6 rows and 4 columns
                   in_tissue array_row array_col   sample_id
                   <integer> <integer> <integer> <character>
AAACAACGAATAGTTC-1         1         0        16    sample01
AAACAAGTATCTCCCA-1         1        50       102    sample01
AAACAATCTACTAGCA-1         0         3        43    sample01
AAACACCAATAACTGC-1         1        59        19    sample01
AAACAGAGCGACTCCT-1         1        14        94    sample01
AAACAGCTTTCAGAAG-1         1        43         9    sample01

Hmmm, check out that in_tissue column! It contains 0/1 because we read in the raw Space Ranger output, which includes all spots including those which don’t actually overlap tissue. Ultimately, we won’t want to analyze the spots that are not on top of tissue, so we’ll have to deal with that.

Having all this image information in the SPE object means we can plot directly. The package ggspavis provides several useful functions for plotting, let’s see the highlights:

# Spots overlaying the H&E, and we can see the slide boundaries as well
ggspavis::plotVisium(spe)

By default this shows the spots, but we can hide them and zoom into the tissue on the slide. This is a helpful way to pop up the H&E for side-by-side comparisons with other plots you might make.

ggspavis::plotVisium(spe, spots = FALSE, zoom = TRUE) 

Let’s take a moment to chat about Wilm’s Tumor - these tumors in the developing kidney are often composed of a couple histologic compartments, blastemal, stromal (with mesenchymal biology), and sometimes epithelial components. In this section, we can see two major components: the bluer indicating densely cellular, ECM-poor regions consistent with blastema, and paler pink regions consistent with stromal tissue.

Let’s actually go ahead and save this plot to a variable; it will be a convenient way for us to quickly pop up the H&E for side-by-side comparison with other plots we’ll make later.

p_he <- ggspavis::plotVisium(spe, spots = FALSE, zoom = TRUE) 

There’s another function called plotCoords which hides the H&E to just show the spots; this plot is not currently very compelling, but once we start coloring by things it will be much more fun!

# no H&E
ggspavis::plotCoords(spe, point_size = 1) 

Quality control

Removing uninformative bins

We only care about spots on top of tissue, so let’s begin by removing the in_tissue = 0 spots. (Note that reading in the filtered_feature_bc_matrix version of Space Ranger output would already be filtered to only 1 in this column, so we’re only taking this step because we read in the raw output).

We can visualize which spots those are, and we’ll do it over the H&E to clearly see the relationship. We’ll use the annotate argument, but plotVisium sees that this column is an integer and forces it to use a continuous color scale; we’ll go ahead and make it a factor for plotting.

# make a factor version of this column to plot with
spe$in_tissue_factor <- as.factor(spe$in_tissue)

ggspavis::plotVisium(
  spe, 
  annotate = "in_tissue_factor", 
  # custom palette so we can see clearly
  pal = c("red", "lightblue")
)

• The tissue overhang on the left isn’t colored at all - indeed, there aren’t spots at those coordinates outside the slide
• Red points are those to filter out - they are uninformative. Worth noting that some of these appear “inside” the tissue.

# keep only spots that are in the tissue
spe <- spe[, spe$in_tissue == 1]
spe

class: SpatialExperiment 
dim: 60664 4210 
metadata(2): resources spatialList
assays(1): counts
rownames(60664): ENSG00000284662 ENSG00000186827 ... ENSG00000277475
  ENSG00000275405
rowData names(3): ID Symbol Type
colnames(4210): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
  TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
colData names(5): in_tissue array_row array_col sample_id
  in_tissue_factor
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):
spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
imgData names(4): sample_id image_id data scaleFactor

Now, we’re down to 3701 spots and we can see this reflected in a plot. We now only see spots with tissue underneath. Onwards!

ggspavis::plotVisium(spe)

Filtering for bin quality

We could do the same kind of QC that we do for single-cell data, which (usually) involves defining global thresholds for e.g. number of detected UMIs or percent mito. However, spatial data doesn’t lend itself as well to this type of filtering since global thresholds may not be suitable for the different tissues. Compared to spatial data, single-cell data more closely (but of course not exactly!) resembles a homogeneous pool of cells.

Let’s go ahead and calculate some QC stats and visualize them:

is_mito <- grepl("(^MT-)|(^mt-)", rowData(spe)$Symbol)
mito_ensembl <- rownames(spe)[is_mito]
spe <- scater::addPerCellQC(spe, subsets=list(mito=mito_ensembl))

ggspavis comes with a helpful QC plotter (makes histograms by default but has a couple more options!)

p1 <- ggspavis::plotObsQC(spe, x_metric = "sum") + ggtitle("sum")
p2 <- ggspavis::plotObsQC(spe, x_metric = "detected") + ggtitle("detected")
p3 <- ggspavis::plotObsQC(spe, x_metric = "subsets_mito_percent") + ggtitle("mito %")

p1 + p2 + p3

Distributions don’t look like there are any major issues with quality here. One reason for this is - these are all spots (aggregates of cells), not cells! Data is sparse, but not quite as sparse.

Let’s plot the same stats on the slide:

p1 <- plotCoords(spe, annotate="sum", point_size = 1) + ggtitle("sum")
p2 <- plotCoords(spe, annotate="detected", point_size = 1) + ggtitle("detected")
p3 <- plotCoords(spe, annotate="subsets_mito_percent", point_size = 1) + ggtitle("mito %")
p1 + p2 + p3 + p_he

We do see a relationship here with tissue type and QC stats that suggests we don’t have globally homogeneous patterns -

• the likely stromal regions look like they have relatively higher mito (although it’s low all around!)
• the likely blastema regions tend to have more detected genes

This tells us that there is local structure in the data, and using global thresholds could be too aggressive/not aggressive enough depending on the local context.

Local filtering with Spot Sweeper

A more spatially-suited approach comes from SpotSweeper which models QC stats to identify local outliers.

We’ll detect local outliers based on these three QC stats. This function will add some colData columns we we can use to find all the local outliers.

spe <- SpotSweeper::localOutliers(spe, metric="sum", direction="lower", log=TRUE) # lower-value outliers should be detected
spe <- SpotSweeper::localOutliers(spe, metric="detected", direction="lower", log=TRUE) # lower-value outliers should be detected
spe <- SpotSweeper::localOutliers(spe, metric="subsets_mito_percent", direction="higher", log=FALSE) # higher-value outliers should be detected

How many of each kind of outlier?

c("sum_outliers", "detected_outliers", "subsets_mito_percent_outliers") |>
  purrr::set_names() |>
  purrr::map(\(x) table(spe[[x]]))

$sum_outliers

FALSE  TRUE 
 4175    35 

$detected_outliers

FALSE  TRUE 
 4169    41 

$subsets_mito_percent_outliers

FALSE  TRUE 
 4190    20 

We’ll use a built-in plotting function from SpotSweeper to plot these results (the ggspavis plotting isn’t quite as nice in this case). Each plot highlights the spots to remove, and spots are colored based on a different stat. Again we’ll include the H&E as a panel for ease of comparison.

p1 <- SpotSweeper::plotQCmetrics(spe, metric = "sum_log", outliers = "sum_outliers") + ggtitle("log sum")
p2 <- SpotSweeper::plotQCmetrics(spe, metric = "detected_log", outliers = "detected_outliers") + ggtitle("log detected")
p3 <- SpotSweeper::plotQCmetrics(spe, metric = "subsets_mito_percent", outliers = "subsets_mito_percent_outliers") + ggtitle("mito %")

p1 + p2 + p3 + p_he

Some of these spots are being flagged by multiple stats, in particular both sum and detected. Let’s go ahead and filter all these cells out:

# combine all outliers into "local_outliers" column to filter on
spe$local_outliers <- as.logical(spe$sum_outliers) |
    as.logical(spe$detected_outliers) |
    as.logical(spe$subsets_mito_percent_outliers)

spe <- spe[, spe$local_outliers == FALSE]

Normalization

With single-cell, a popular normalization approach we adopt involves deconvolution - we do a quick clustering so we can get groups of cells that will share a size factor, for faster processing. But this isn’t suitable for spot-based spatial data since our units aren’t cells, but potentially heterogeneous groups of cells (whole or partial).

We’ll therefore use going to use library-size normalization here to correct for sequencing depth but preserve the other heterogeneous biological signal. Bear in mind the caveat that in some cases library size itself is [also confounded with spatial structure]((https://link.springer.com/article/10.1186/s13059-024-03241-7); so, keep your eyes peeled for normalization developments, since this space is evolving rapidly!

spe <- scuttle::computeLibraryFactors(spe)
spe <- scuttle::logNormCounts(spe)
spe # now we have a logcounts assay

class: SpatialExperiment 
dim: 60664 4149 
metadata(2): resources spatialList
assays(2): counts logcounts
rownames(60664): ENSG00000284662 ENSG00000186827 ... ENSG00000277475
  ENSG00000275405
rowData names(4): ID Symbol Type subsets_mito
colnames(4149): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
  TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
colData names(21): in_tissue array_row ... local_outliers sizeFactor
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):
spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
imgData names(4): sample_id image_id data scaleFactor

Dimension reduction

For this step, we’re going to get a reduced dimension representation using just the lens of gene expression, “ignoring” spatial information. (Later, we’ll see complementary approaches to some of this which leverages spatial information).

This all tells us we have to be careful when interpreting analyses on spatial data: “to what extent was spatial information taken into account” is something you’ll want to understand.

Let’s do HVG selection (again, this does not consider spatial patterns! we’ll learn about spatially variable genes later), PCA, UMAP, and clustering (once again, this won’t be spatial regions! just spot similarity based on expression - that’s not the same as cells, either!) So, we’ll use a SCE-style approach:

num_genes <- 2000
gene_variance <- scran::modelGeneVar(spe)
hv_genes <- scran::getTopHVGs(gene_variance, n = num_genes)

spe <- scater::runPCA(spe, subset_row = hv_genes)
spe <- scater::runUMAP(spe, dimred = "PCA")

nn_clusters <- scran::clusterCells(
  spe, # object to perform clustering on
  use.dimred = "PCA", # perform clustering on the PCA matrix
  BLUSPARAM = bluster::NNGraphParam(
    k = 20,
    type = "jaccard",
    cluster.fun = "louvain"
  )
)
spe$nn_cluster <- nn_clusters

We’ll plot the UMAP colored by clusters, because we love UMAPs:

scater::plotUMAP(spe, color_by = "nn_cluster") +
  # use distinct palette so we can match up with a spatial view next
  ggokabeito::scale_color_okabe_ito()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

How do these clusters come out when using the spatial layout? We have coordinates, after all - let’s use them!

p1 <- ggspavis::plotCoords(spe, annotate = "nn_cluster", point_size = 1) +
  # matching palette to UMAP
  ggokabeito::scale_color_okabe_ito()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

p1 + p_he

A lot of the clusters are intermingled in space, showing that this analysis is really not spatially-aware! But, we do see some interesting patterns - check out clusters 1, 5, & 7 in particular - they look stroma-y; also, cluster 8 ooks like a defined region in the H&E, but I’m not a pathologist so let’s not overinterpret here. When expression and tissue structure are correlated, expression-only clusters can pick up on certain spatial regions, but this is not the same as finding “true” spatial domains.

Exploring spot-based spatial data

Let’s explore this data a little more, using some marker genes. We don’t have specific cells as units here, but we can use marker gene expression to build some descriptive intuition for the data.

We’ve defined some marker genes here for visualization that correspond to expected Wilm’s Tumor biology:

marker_genes <- c(
  "WT1 (kidney development, maglinant in WT)" = "ENSG00000184937",
  "SIX1 (blastema, malignant)" = "ENSG00000126778",
  "MYCN (malignant)" = "ENSG00000134323",
  "COL1A1 (fibroblast/mesenchymal)" = "ENSG00000108821",
  "TAGLN (myofibroblast/mesenchymal)" = "ENSG00000149591"
)

Let’s plot the expression of these genes, using both UMAP and the spatial layouts to compare:

cluster_umap <- scater::plotUMAP(spe, color_by = "nn_cluster", point_size = 0.5, point_alpha = 0.5) +
  ggtitle("expression-based clusters") + 
  coord_equal() +
  ggokabeito::scale_color_okabe_ito() 

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

marker_umaps <- marker_genes |>
  purrr::imap(
    \(ensembl, gene_name) {
      scater::plotUMAP(spe, color_by = ensembl, point_size = 0.5, point_alpha = 0.5) + ggtitle(gene_name) + coord_equal()
    }
  ) |>
  # cluster umap and then the marker gene umaps
  append(values = cluster_umap, after = 0)

patchwork::wrap_plots(marker_umaps, nrow = 2)

cluster_coords <- ggspavis::plotCoords(spe, annotate = "nn_cluster", point_size = 0.75) +
  ggtitle("expression-based clusters") + 
  ggokabeito::scale_color_okabe_ito() 

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

marker_coords <- marker_genes |>
  purrr::imap(
    \(ensembl, gene_name) {
      # default assay is counts, override!
      ggspavis::plotCoords(spe, annotate = ensembl, point_size = 0.75, assay = "logcounts") +
        scale_color_distiller(palette = "Blues", direction = 1) + 
        ggtitle(gene_name)
    }
  ) |>
  # cluster cooords and then the marker gene coords
  append(values = cluster_coords, after = 0) 

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

p_he + patchwork::wrap_plots(marker_coords, nrow = 2) + plot_layout(widths = c(1,2))

Based on these observations, it seems like there could be some spatial patterns/variability here; interpret. We’ll talk more rigorously about quantifying spatially-variable genes later in the workshop, but for now, we can take this opportunity to to ask some descriptive questions about these genes’ expression across space.

Spatial descriptive statistics

Although spatial analysis is an emerging field in transcriptomics, it’s not new to science! Folks in geography, epidemiology, ecology, and more have been thinking about spatial data for a long time, so we can draw from existing quantitative approaches.

Spot-based spatial data is analogous to “lattice-structure” spatial data, in contrast to high-res spatial technologies that are “point-pattern” data; see the pasta paper.

The R package spdep already exists to work with lattice-structure spatial data, and helpfully folks have already started to put this (and other spatial platforms) to work for transcriptomics data. The Voyager/SpatialFeatureExperiment: https://pachterlab.github.io/voyager/index.html framework is complementary to SpatialExperiment and can help us use these tools. The SFE object is an expanded SPE object with more slots to facilitate spatial analysis with existing tools (we won’t get into these weeds much though, suffice to say it’s a super rich environment).

When working with spot data, we can turn it into a graph and ask descriptive (not modeling!) questions about spatial structure directly. Let’s go!

First, we’ll need to get it into the SFE framework and convert the SPE to an SFE object.

# TODO: I have to remove the image before converting, errors otherwise. Can this be avoided?
# Error: [ext] invalid extent
spe_noimg <- spe
imgData(spe_noimg) <- imgData(spe_noimg)[0, , drop=FALSE]
sfe <- toSpatialFeatureExperiment(spe_noimg)
#sfe <- mirror(sfe, direction = "vertical") # use mirror so plots between ggspavis and this one match. 
# OR, let it be a mirror to serve as a teachable moment about space!

One way to make a graph is using k-nearest neighbors, but there’s other ways (all of these choices matter to some degree, and we should discuss that!): Now we’ll make a graph using:

• spot centroids as nodes (don’t skip that this is a choice!)
• knn approach with k = 10, as a starting point
• no weighting

colGraph(sfe, "knn10") <- SpatialFeatureExperiment::findSpatialNeighbors(
  sfe, 
  type = "centroids", 
  method = "knearneigh", 
  k = 10
)

This will get stored in colGraph slot in the SFE object.

We can visualize it - check out those edges and nodes, it’s a graph!

Voyager::plotColGraph(
  sfe, 
  colGraphName = "knn10", 
  colGeometryName = "centroids", 
  segment_size = 0.1,
  geometry_size = 0.3
) +
  theme_void() # we just want the graph

Hey notice anything interesting about this image? It’s a mirror of how ggspavis plotted it! Is that a problem? Well, no - in space we can mirror and rotate and such, but the spatial relationships are preserved. This is a big hint that specific coordinates may not always be the same, but the distances should be.

It seems potentially cute to zoom in to prove to you that k = 10:

Voyager::plotColGraph(
  sfe, 
  colGraphName = "knn10", 
  colGeometryName = "centroids", 
  segment_size = 0.1,
  geometry_size = 5
) +
  scale_y_continuous(limits = c(11500, 12000)) + 
  scale_x_continuous(limits = c(4000,4300)) +
  theme_void()

Alright, now that we have a graph, let’s see an example of using it for exploratory analysis. A common descriptive statistic used for lattice spatial data is “Moran’s I” which quantifies spatial autocorrelation. It asks, is the spatial distribution of a variable of interest dispersed, random, or clustered? We’ll use it to look at gene expression of some of our marker genes of interest. In this case, Moran’s I tells us: do neighboring spots (hint: this will depend on your definition of a neighborhood!_) have very similar, very dissimilar, or totally unrelated gene expression, where:

• close to 1 values: Nearby spots tend to have similar expression
• ~0: looks random in space
• close to -1 values: Nearby spots tend to have different expression (less commonly of interest in this sort of biology)

We’ll calculate this with Voyager on our SFE object using the graph we just built.

sfe <- Voyager::runUnivariate(
  sfe, 
  type = "moran", 
  features = marker_genes, 
  colGraphName = "knn10"
)

This stores results in rowData:

rowData(sfe)[marker_genes, ]

DataFrame with 5 rows and 6 columns
                             ID      Symbol            Type subsets_mito
                    <character> <character>     <character>    <logical>
ENSG00000184937 ENSG00000184937         WT1 Gene Expression        FALSE
ENSG00000126778 ENSG00000126778        SIX1 Gene Expression        FALSE
ENSG00000134323 ENSG00000134323        MYCN Gene Expression        FALSE
ENSG00000108821 ENSG00000108821      COL1A1 Gene Expression        FALSE
ENSG00000149591 ENSG00000149591       TAGLN Gene Expression        FALSE
                moran_sample01 K_sample01
                     <numeric>  <numeric>
ENSG00000184937       0.400004    5.36664
ENSG00000126778       0.146762    2.52327
ENSG00000134323       0.116602    2.06943
ENSG00000108821       0.598150    2.90257
ENSG00000149591       0.768734    4.80837

This makes some sense given what we when plotting gene expression:

• SIX1 and MYCN were just splotched all around randomly in much of the slide
• COL1A1 and TAGLN were more concentrated in certain areas
• WT1 was somewhere in the middle - generally expressed by less variation from spot-to-spot in the plots earlier, per vibes

This global statistic gave us an overall picture, but we can calculate a variant of this statistic also called Local Moran’s I which is, you guessed it, a local version. Rather than giving a single value for the dataset, we’ll get a per-spot value (or rather, per node in the graph) that tells us “how much does this spot agree with its immediate neighborhood”?

• high: spot is similar to neighbors and all are high OR all are low
• ~0: spatially neutral, no real relationship to neighbors
• low: spot is different from neighbors where high surrounded by low or vice versa

Let’s go:

sfe <- Voyager::runUnivariate(
  sfe, 
  type = "localmoran", 
  features = marker_genes, 
  colGraphName = "knn10", 
  name = "localmoran_knn10" # save it as this name
)

# stores it here, not so pretty...
localResults(sfe)$localmoran_knn10

DataFrame with 4149 rows and 5 columns
                                             ENSG00000184937
                                                <data.frame>
AAACAACGAATAGTTC-1  1.667219796:-3.16454e-04:1.30970e-01:...
AAACAAGTATCTCCCA-1  0.007633311:-1.74338e-07:7.21757e-05:...
AAACACCAATAACTGC-1 -0.000607848:-1.18586e-09:4.90946e-07:...
AAACAGAGCGACTCCT-1  0.051163178:-2.43975e-04:1.00981e-01:...
AAACAGCTTTCAGAAG-1  0.078470169:-3.93378e-05:1.62852e-02:...
...                                                      ...
TTGTTTCACATCCAGG-1    -0.1995113:-5.79632e-04:0.23982823:...
TTGTTTCATTAGTCTA-1    -0.0648147:-8.20217e-05:0.03395418:...
TTGTTTCCATACAACT-1     0.0496851:-5.42162e-05:0.02244428:...
TTGTTTGTATTACACG-1     0.0765047:-1.90444e-05:0.00788423:...
TTGTTTGTGTAAATTC-1     2.1740629:-3.44664e-04:0.14264156:...
                                            ENSG00000126778
                                               <data.frame>
AAACAACGAATAGTTC-1 -2.15283539:-3.16619e-04:0.131038522:...
AAACAAGTATCTCCCA-1  0.12587206:-9.83964e-05:0.040732072:...
AAACACCAATAACTGC-1 -0.00824221:-7.23225e-07:0.000299415:...
AAACAGAGCGACTCCT-1 -0.52228726:-3.17655e-04:0.131467282:...
AAACAGCTTTCAGAAG-1 -0.23168190:-1.34973e-04:0.055871160:...
...                                                     ...
TTGTTTCACATCCAGG-1   -0.0939208:-1.35766e-05:0.00562064:...
TTGTTTCATTAGTCTA-1   -0.0682406:-5.91852e-04:0.24488153:...
TTGTTTCCATACAACT-1    0.0451762:-6.60282e-05:0.02733382:...
TTGTTTGTATTACACG-1    0.1783573:-7.46413e-05:0.03089916:...
TTGTTTGTGTAAATTC-1   -0.0634922:-6.58503e-06:0.00272618:...
                                          ENSG00000134323
                                             <data.frame>
AAACAACGAATAGTTC-1 -0.8411838:-1.16198e-04:0.04810041:...
AAACAAGTATCTCCCA-1 -0.0491596:-5.23386e-04:0.21656837:...
AAACACCAATAACTGC-1 -0.0717147:-4.11319e-05:0.01702787:...
AAACAGAGCGACTCCT-1 -0.0755661:-1.11880e-05:0.00463178:...
AAACAGCTTTCAGAAG-1 -0.0174028:-2.23178e-05:0.00923935:...
...                                                   ...
TTGTTTCACATCCAGG-1 0.0111875:-1.11170e-06:0.000460242:...
TTGTTTCATTAGTCTA-1 0.2150468:-7.14312e-05:0.029570378:...
TTGTTTCCATACAACT-1 0.2422880:-2.46435e-04:0.101998874:...
TTGTTTGTATTACACG-1 0.0572942:-1.06836e-04:0.044225137:...
TTGTTTGTGTAAATTC-1 0.3536584:-5.23386e-04:0.216568367:...
                                             ENSG00000108821
                                                <data.frame>
AAACAACGAATAGTTC-1      3.097516:-1.00038e-03:0.41374280:...
AAACAAGTATCTCCCA-1     -0.049879:-9.91574e-06:0.00410507:...
AAACACCAATAACTGC-1      0.159551:-1.96079e-05:0.00811751:...
AAACAGAGCGACTCCT-1     -0.140422:-3.45277e-05:0.01429394:...
AAACAGCTTTCAGAAG-1      0.321191:-9.37702e-05:0.03881716:...
...                                                      ...
TTGTTTCACATCCAGG-1 -0.159006137:-2.35465e-05:9.74803e-03:...
TTGTTTCATTAGTCTA-1  0.104516002:-7.92831e-05:3.28206e-02:...
TTGTTTCCATACAACT-1  0.000407591:-2.20039e-09:9.10962e-07:...
TTGTTTGTATTACACG-1  1.949902019:-1.46119e-03:6.04047e-01:...
TTGTTTGTGTAAATTC-1  0.848563795:-2.20292e-04:9.11805e-02:...
                                           ENSG00000149591
                                              <data.frame>
AAACAACGAATAGTTC-1  5.2638003:-1.08466e-03:0.448560970:...
AAACAAGTATCTCCCA-1 -0.0363059:-1.48510e-06:0.000614829:...
AAACACCAATAACTGC-1  0.2883653:-1.40370e-04:0.058104774:...
AAACAGAGCGACTCCT-1 -0.0448511:-1.64269e-05:0.006800604:...
AAACAGCTTTCAGAAG-1  0.4025505:-1.40370e-04:0.058104774:...
...                                                    ...
TTGTTTCACATCCAGG-1   0.3590522:-1.40370e-04:0.05810477:...
TTGTTTCATTAGTCTA-1   0.0369287:-4.97746e-06:0.00206066:...
TTGTTTCCATACAACT-1  -0.0548757:-8.82533e-06:0.00365365:...
TTGTTTGTATTACACG-1   0.1463867:-1.31907e-04:0.05460239:...
TTGTTTGTGTAAATTC-1   2.6481197:-8.84656e-04:0.36592308:...

# we can plot with voyager, since it's in a special slot
# or we'd have to pull it out to store back in our SPE
Voyager::plotLocalResult(
  sfe,
  name = "localmoran_knn10",
  features = marker_genes,
  colGeometryName = "centroids",
  divergent = TRUE, # palette
  diverge_center = 0 # palette 
  ) 

For easier comparison, we can also use the standardized (z-score) values by specifying the specific attribute (the default was "Ii", the actual stat).

Voyager::plotLocalResult(
  sfe,
  name = "localmoran_knn10",
  features = marker_genes,
  attribute = "Z.Ii", # plot the normalized moran's i
  colGeometryName = "centroids",
  divergent = TRUE, 
  diverge_center = 0 
  ) 

Scale matters!

We did a graph with 10 nearest neighbors. The number of neighbors is going to matter here!

Pick one gene (provide the ensembl id) and see how interpretation differs at different k’s. Try out a really high (k>=50) and really low number (k<=3) below to calculate and plot Local Moran’s I again - the real value, not scaled. What do you see?

gene <- marker_genes[1] # WT1, for example
colGraph(sfe, "knn100") <- findSpatialNeighbors(sfe, type="centroids", method = "knearneigh", k=100)
colGraph(sfe, "knn2") <- findSpatialNeighbors(sfe, type="centroids", method = "knearneigh", k=2)

sfe <- Voyager::runUnivariate(
  sfe, 
  type="localmoran", 
  features = gene, 
  colGraphName="knn100", 
  name = "localmoran_knn100" 
)
sfe <- Voyager::runUnivariate(
  sfe, 
  type="localmoran", 
  features = gene, 
  colGraphName="knn2", 
  name = "localmoran_knn2" 
)

p1 <- Voyager::plotLocalResult(
  sfe,
  name = "localmoran_knn100",
  features = gene,
  colGeometryName = "centroids",
  divergent = TRUE,
  diverge_center = 0
) 


p2 <- Voyager::plotLocalResult(
  sfe,
  name = "localmoran_knn2",
  features = gene,
  colGeometryName = "centroids",
  divergent = TRUE,
  diverge_center = 0
) 

p1 + p2

When k is very high, there’s a lot of neighbors being consider. Expression will have to be similar cross them all to get a high value for this stat. But, when k is very low, there’s very few neighbors! Much finer-grained picture of similarity between spots.

Session Info

sessionInfo()

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices datasets  utils     methods  
[8] base     

other attached packages:
 [1] Voyager_1.12.0                  SpatialFeatureExperiment_1.12.1
 [3] patchwork_1.3.2                 ggspavis_1.16.0                
 [5] ggplot2_4.0.1                   VisiumIO_1.6.3                 
 [7] TENxIO_1.12.1                   SpatialExperiment_1.20.0       
 [9] SingleCellExperiment_1.32.0     SummarizedExperiment_1.40.0    
[11] Biobase_2.70.0                  GenomicRanges_1.62.1           
[13] Seqinfo_1.0.0                   IRanges_2.44.0                 
[15] S4Vectors_0.48.0                BiocGenerics_0.56.0            
[17] generics_0.1.4                  MatrixGenerics_1.22.0          
[19] matrixStats_1.5.0              

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.23          splines_4.5.2            
  [3] BiocIO_1.20.0             bitops_1.0-9             
  [5] tibble_3.3.1              R.oo_1.27.1              
  [7] lifecycle_1.0.5           sf_1.0-24                
  [9] edgeR_4.8.2               rprojroot_2.1.1          
 [11] lattice_0.22-7            MASS_7.3-65              
 [13] magrittr_2.0.4            limma_3.66.0             
 [15] sass_0.4.10               rmarkdown_2.30           
 [17] jquerylib_0.1.4           yaml_2.3.12              
 [19] metapod_1.18.0            otel_0.2.0               
 [21] ggside_0.4.1              sp_2.2-0                 
 [23] DBI_1.2.3                 RColorBrewer_1.1-3       
 [25] multcomp_1.4-29           abind_1.4-8              
 [27] spatialreg_1.4-2          purrr_1.2.1              
 [29] R.utils_2.13.0            RCurl_1.98-1.17          
 [31] TH.data_1.1-5             sandwich_3.1-1           
 [33] ggrepel_0.9.6             irlba_2.3.5.1            
 [35] terra_1.8-93              units_1.0-0              
 [37] RSpectra_0.16-2           dqrng_0.4.1              
 [39] DelayedMatrixStats_1.32.0 codetools_0.2-20         
 [41] DropletUtils_1.30.0       DelayedArray_0.36.0      
 [43] scuttle_1.20.0            tidyselect_1.2.1         
 [45] memuse_4.2-3              farver_2.1.2             
 [47] ScaledMatrix_1.18.0       viridis_0.6.5            
 [49] jsonlite_2.0.0            escheR_1.10.0            
 [51] BiocNeighbors_2.4.0       e1071_1.7-17             
 [53] survival_3.8-3            scater_1.38.0            
 [55] tools_4.5.2               ggnewscale_0.5.2         
 [57] Rcpp_1.1.1                glue_1.8.0               
 [59] gridExtra_2.3             SparseArray_1.10.8       
 [61] BiocBaseUtils_1.12.0      xfun_0.56                
 [63] here_1.0.2                EBImage_4.52.0           
 [65] dplyr_1.1.4               HDF5Array_1.38.0         
 [67] withr_3.0.2               BiocManager_1.30.27      
 [69] fastmap_1.2.0             boot_1.3-32              
 [71] rhdf5filters_1.22.0       bluster_1.20.0           
 [73] spData_2.3.4              digest_0.6.39            
 [75] rsvd_1.0.5                R6_2.6.1                 
 [77] wk_0.9.5                  LearnBayes_2.15.1        
 [79] jpeg_0.1-11               R.methodsS3_1.8.2        
 [81] h5mread_1.2.1             renv_1.1.6               
 [83] data.table_1.18.0         class_7.3-23             
 [85] htmlwidgets_1.6.4         S4Arrays_1.10.1          
 [87] spdep_1.4-1               uwot_0.2.4               
 [89] pkgconfig_2.0.3           scico_1.5.0              
 [91] gtable_0.3.6              S7_0.2.1                 
 [93] XVector_0.50.0            htmltools_0.5.9          
 [95] fftwtools_0.9-11          scales_1.4.0             
 [97] png_0.1-8                 scran_1.38.0             
 [99] SpotSweeper_1.6.0         knitr_1.51               
[101] tzdb_0.5.0                rjson_0.2.23             
[103] coda_0.19-4.1             nlme_3.1-168             
[105] proxy_0.4-29              cachem_1.1.0             
[107] zoo_1.8-15                rhdf5_2.54.1             
[109] KernSmooth_2.23-26        parallel_4.5.2           
[111] vipor_0.4.7               s2_1.1.9                 
[113] pillar_1.11.1             grid_4.5.2               
[115] vctrs_0.7.0               BiocSingular_1.26.1      
[117] beachmat_2.26.0           sfheaders_0.4.5          
[119] cluster_2.1.8.1           beeswarm_0.4.0           
[121] evaluate_1.0.5            readr_2.1.6              
[123] zeallot_0.2.0             magick_2.9.0             
[125] mvtnorm_1.3-3             cli_3.6.5                
[127] locfit_1.5-9.12           compiler_4.5.2           
[129] rlang_1.1.7               labeling_0.4.3           
[131] classInt_0.4-11           spatialEco_2.0-3         
[133] ggbeeswarm_0.7.3          viridisLite_0.4.2        
[135] deldir_2.0-4              BiocParallel_1.44.0      
[137] tiff_0.1-12               Matrix_1.7-4             
[139] hms_1.1.4                 sparseMatrixStats_1.22.0 
[141] Rhdf5lib_1.32.0           statmod_1.5.1            
[143] igraph_2.2.1              bslib_0.9.0              
[145] ggokabeito_0.1.0